Collective behavior of heterogeneous neural networks 
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We investigate a network of integrate-and-fire neurons characterized by a distribution of spik- 
ing frequencies. Upon increasing the coupling strength, the model exhibits a transition from an 
asynchronous regime to a nontrivial collective behavior. At variance with the Kuramoto model, (i) 
the macroscopic dynamics is irregular even in the thermodynamic limit, and (ii) the microscopic 
(single-neuron) evolution is linearly stable. 
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The investigation of networks of oscillators can pro- 
vide new insights on the basic mechanisms which under- 
lie brain functioning. In particular, the spontaneous on- 
set of a collective dynamics is an intriguing phenomenon 
that can contribute to information transmission across 
different brain areas. Given the large number N of neu- 
rons (oscillators) present in a real brain, it is tempting to 
adopt a statistical-mechanics point of view and thereby 
investigate the behavior for N — > oo (the so-called ther- 
modynamic limit). Two different setups are typically in- 
voked [1]: (i) sparse networks, characterized by a fixed 
number of synaptic connections; (ii) massively connected 
networks, where the number of connections is propor- 
tional to N. In the former case, the local field seen by 
the single neurons naturally fluctuates even for N — > oo 
(being the sum of a fixed finite number of different in- 
put signals) , consistently with the experimental evidence 
of an irregular background activity in the cerebral cor- 
tex . The latter setup has the advantage of being, at 
least in principle, amenable to an exact mean-field treat- 
ment, although microscopic fluctuations survive only if 
inhibition and excitation balance each other Q. 

In this Letter, we numerically show that an irregu- 
lar microscopic and macroscopic dynamics can generi- 
cally arise even in an inhibitory, globally coupled net- 
work. More precisely, we consider a heterogeneous net- 
work of pulse-coupled integrate-and-fire (IF) neurons |4| , 
each characterized by a different bare spiking frequency. 
This setup is similar to that of the Kuramoto model 
(KM) [51, where each single oscillator is identified by a 
phase variable 4>. The analogy is so tight that it has even 
been shown that the pulse-coupling mechanism charac- 
terizing IF neurons reduces, in the weak coupling limit, 
to that of the KM, the only difference being that the 
coupling function is not purely sinusoidal Q . It is there- 
fore quite important to clarify to what extent a network 
of IF neurons reproduces the KM scenario for stronger 
coupling strengths, especially by recalling that the KM 
is often invoked while testing new ideas on the control 
of synchronization within neural contexts 0] ■ Finally, in 
order to make the model closer to a realistic setup, we 



include delay to account for the finite propagation time 
of the electric pulses. 

Our strategy consists in studying the macroscopic col- 
lective dynamics in the large JV-limit, for different values 
of the coupling strength g. In the KM, it is known that 
for a weak enough coupling, the single oscillators rotate 
independently of each other. On the other hand, above 
a critical value, a subset of them mutually synchronize, 
as signalled by a non zero value of the order parameter 
p = |(e l ^)| (the angular brackets denote an average over 
the rotators). In this Letter we show that IF neurons 
give rise to a similar but substantially different scenario. 
First of all, the (second) maximum Lyapunov exponent 
is always negative implying that the evolution must 
eventually converge to a periodic orbit. On the other 
hand, the study of relatively small networks shows that 
the time needed to approach a periodic orbit is expo- 
nentially long with the system size, implying that the 
"transient" extends over increasingly longer time scales. 
In other words, this is an instance of stable chaos 0, a 
phenomenon already detected in networks of pulse cou- 
pled oscillators without delay [l(|, although its onset in 
systems with delayed coupling is controversial 
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A second difference is that, at variance with the KM, 
the coupling contributes also to slowing down the spiking 
activity of the single neurons (a somehow similar mecha- 
nism operates in ensemble of cold atoms [Hj]) and drives 
a subset of neurons below the firing threshold - a phe- 
nomenon reminiscent of oscillator-death However, 
the most striking difference concerns the above-threshold 
regime, as the overall neural activity is not simply pe- 
riodically modulated, but exhibits irregular, seemingly 
chaotic, oscillations (still in the presence of a negative 
"microscopic" second Lyapunov exponent). Nothing of 
this type has been observed in the corresponding setup 
of a KM with delayed coupling 15j . 

The evolution equations for the N membrane poten- 
tials Vi write, 
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where all variables are expressed in adimensional units. 
When Vi reaches the threshold u, = 1, it is instan- 
taneously reset to the value Vi = 0, while a spike is 
emitted (and received with a delay t c i). The network 
is assumed to be heterogeneous, in that different neu- 
rons are exposed to different suprathreshold currents af, 
(^o)i = l/ln[(aj/(ai — 1)] is the bare spiking frequency. 
Sij denotes the connectivity matrix and the sum in 
Eq. ([l]) runs over the spikes received by the neuron i. 
Finally, the coupling strength g is our control parameter: 
the negative sign in front of last term in the r.h.s. means 
that we assume inhibitory coupling. Notice also that the 
same last term does not only couple the oscillators but 
modifies also their frequency. 

All the simulations reported in this Letter refer to a 
globally coupled network, i.e. SVz = 1 for any i, I, but we 
have verified that the introduction of additional disorder 
(by randomly removing a fixed fraction of connections) 
does not substantially modify the overall scenario. The 
delay is set everywhere equal to t<j = 0.1, while the cur- 
rents cii are randomly and uniformly distributed in the 
interval [1.2,2.8]. These parameter values are consistent 
with those selected in Ref . [Tljj , where they have been cho- 
sen on the basis of biological motivations. In our case, 
the Kuramoto order parameter p cannot be used to char- 
acterize the onset of a collective dynamics. In fact, the 
inhibitory coupling may drive the potential below the 
reset value, thus making the transformation of the Vi po- 
tential into a phase-like variable ill-defined. The difficulty 
can be overcome by coarse graining the spiking activity. 
We do so by dressing each spike with a finite width and 
thereby construct a smooth effective field E. If we as- 
sume the pulse shape, pit) := a 2 iexp(— at) (i > 0), the 
corresponding field E can be generated by integrating the 
equation, 
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This procedure is often used to determine the field ac- 
tually seen by the single neurons [l(|; here it is just 
a strategy to construct a meaningful order parameter, 
that is defined as the standard deviation a of E (er 2 = 
(E 2 ) t — (E)t, where (-)t denotes a time average). We 
choose a = 20, a value that corresponds to sufficiently 
broad pulses to get rid of the statistical fluctuations, but 
not so large as to wash out the time evolution. As long 
as the asymptotic regime is an asynchronous state char- 
acterized by a constant activity, a is zero in the infinite 
N limit, while any form of collective dynamics gives rise 
to a nonzero a. This is precisely what is seen in Fig. [T] 
where a is plotted versus the coupling strength g for dif- 
ferent network sizes. Below g c ss 0.5, a is quite small and 
appears to decrease as 1/VN with the system size (see 
the left inset), indicating that the deviation from zero 
is a finite-size effect. Above g c , a starts to grow and is 



independent of the system size, suggesting the onset of 
some form of synchronization (the right inset contains an 
instance of the field evolution for g = 5). Superficially, 
this scenario is reminiscent of the synchronization tran- 
sition observed in the KM. In the following we show that 
there arc several conceptually relevant differences. The 




FIG. 1: (Color online) Standard deviation, a, of the effective 
field £ as a function of the coupling strength g for iV=5,750, 
(red) squares, jV=ll,500, (black) triangles, and 7V=46,000, 
(green) circles. The upper inset contains the rescaled stan- 
dard deviation. The lower inset contains an instance of the 
time evolution of E for g = 5 and N =5,750. 



first difference concerns the microscopic (single neuron) 
behavior. The maximum Lyapunov exponent, A, of the 
Poincare map (to get rid of the first zero Lyapunov ex- 
ponent) is negative both below and above the transition 
and does not depend on N for large system sizes 
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Altogether, the stable microscopic dynamics observed in 
this setup contrasts with the weakly unstable dynamics 
observed in the KM, where the maximum Lyapunov ex- 
ponent is positive, though scales as 1/N [18|. On the 
other hand, the transient time T r needed for a generic 
trajectory to converge to some periodic orbit grows expo- 
nentially with N . This is illustrated in Fig. [5J where the 
average T r (over more than 100 realizations of the dis- 
order) is plotted for different coupling strengths. There, 
one can also appreciate that the exponential growth rate 
decreases systematically with increasing g. Therefore, for 
large N, the relevant dynamical regime is represented by 
the "transient" dynamics, rather than by the periodic or- 
bit approached over astronomical time scales. This stable 
chaos scenario was first observed in the absence of de- 
lay [13] for networks of identical oscillators, when disor- 
der in the connectivity matrix is included. Its occurrence 
in the presence of delay is somehow controversial. It ap- 
pears that the length of the transient depends crucially 
on the balance between the amplitude of the effective 
disorder and the stability of clustered states. Whenever 
local fluctations decrease with N, the transient length 
does not only stops growing exponentially, but even de- 
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FIG. 3: Return maps for the maxima of the effective field E 
when g = 5 for AT=11,500 (a) and 7V=46,000 (b). 



FIG. 2: (Color online) Average transient length T r as a func- 
tion of N for g = 0.3, 0.7, 1.3, 3, and 5 (diamonds, circles, tri- 
angles, squares, and plusses, respectively). The (red) dashed 
lines are exponential fits. 



creases, since generic trajectories rapidly approach one of 
the clustered state s fltl . At variance with the previously 
considered setups [ll|, [l2[, the disorder induced by the 
heterogeneity of the currents, survives in the thermody- 
namic limit. Accordingly, the exponential growth of the 
transient length is expected to persist for arbitrarily large 
N even in the absence of disorder in the connections. In 
fact, we find no evidence of a convergence towards more 
coherent states. 

The standard deviation a allows identifying the very 
existence of collective fluctuations, but does not tell us 
anything about their dynamical character. Up to g ss 2, 
simulations performed for increasing TV suggest that the 
field E behaves periodically in the thermodynamic limit. 
On the other hand, the right inset in Fig. [TJ which refers 
to g = 5, reveals a rather irregular behavior still for 
N ^ 46, 000. A more accurate analysis is however nec- 
essary before making any claim. As a first test, we con- 
struct a return map by plotting the (n + l)-st maximum 
EM(n + 1) of the field versus the previous one. In Fig. [3J 
we see that the points in the Poincare section fill a broad 
and almost the same area for both N =11,500 and 46,000. 
Such features consistently indicate that the collective dy- 
namics is characterized by complex oscillations. 

Next we characterize the collective motion by com- 
puting the Fourier power spectrum S(u) of the field E. 
The spectra reported in Fig. [4] reveal several broad peaks 
whose width does not appear to decrease for increasing 
N. This confirms that the irregularity of the collective 
dynamics persists in networks of arbitrary size and there- 
fore differs from the periodic oscillations reported, e.g. in 
20|,[2l|. 




Rcfs. 

In order to shed further light on the system evolu- 
tion, we have analysed the single-neuron behavior too. 
In Fig. [5J the spiking frequency v (defined as the inverse 



FIG. 4: (Color online) Power spectrum S(y) of the effective 
field for N = 5, 750, (dotted line) 11,500 (dashed line) and 
46,000 (solid line) . The spectra have been obtained by Fourier 
transforming a signal of temporal length w 49, 300. The ar- 
rows point to the frequencies identified in Fig. [5] 



of the average inter-spike interval - ISI) of the single neu- 
rons is plotted versus the bare frequency vq £ [0.558, 2.26] 
(again for g = 5). The effective frequency is systemati- 
cally smaller than vq] this is because the inhibitory cou- 
pling lessens the neural activity. In fact, inhibition is so 
strong, as to bring the least active neurons below thresh- 
old: neurons with i>q <sa 1.56 do not fire at all, and thus 
do not actively contribute to the network dynamics: they 
are just slaved by the other degrees of freedom. 

The appearance of plateaux (the widest ones corre- 
sponding to harmonics of the frequency v — 0.23) reveals 
that neurons with similar bare frequencies lock together, 
as it is naturally expected for periodically forced oscilla- 
tors (see the phenomenon of Arnold tongues). However, 
in this case, the forcing field is not periodic: the shaded 
region around the curve v{vq) highlights the fluctuations 
of the ISI (its vertical width is equal to three standard 
deviations of v). We have verified that such a width does 
not vanish upon increasing N , while the neurons within 
the same plateau are frequency- but not phase-locked. It 
is tempting to trace back the irregular collective motion 
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FIG. 5: (Color online) Average spiking frequency of the os- 
cillators as a function of the bare frequency (ranging in the 
interval [0.558, 2.26]), for g = 5 and N =46,000 neurons (red 
solid line). The thin horizontal lines are located at multiples 
of v — 0.23. The shaded area identifies the region covered by 
frequency fluctuations (see the text for further details). 



to the presence of neurons that are nearly at threshold, 
whose activity is quite sporadic. However, we have veri- 
fied that the overall evolution is almost unchanged when 
such neurons (and those which do not spike at all) are 
removed from the outset. 

All of our numerical simulations confirm that the irreg- 
ularity of the collective dynamics persists for N — » oo. 
It is important to realize that this scenario is a priori 
legitimate, since the dynamics is ruled, in the thermo- 
dynamic limit, by a suitable nonlinear functional equa- 
tion. In this case, the relevant object is the probabil- 
ity density P{v, vq, t) for the membrane potential of the 
neurons, whose bare spiking frequency lies in the interval 
[i/q, i/Q + diso] to belong to the interval [v, v + dv] at time t. 
As functional equations involve infinitely many degrees of 
freedom, one can, in principle, expect an arbitrary degree 
of dynamical complexity. In models such as the networks 
considered in the corresponding probability density 
is a Gaussian and it is therefore described by just two 
variables. As a result, in that context one cannot ob- 
serve anything more complex than periodic oscillations. 
In the standard KM it has been proved that not even 
periodic oscillations can arise; a periodic collective mo- 
tion can be observed only by invoking a more complicate 



nonlinear dependence on the order parameter 22j. On 
the other hand, globally coupled logistic maps exhibit an 
infinite dimensional dynamics [23j |. The problem of de- 
termining the active modes in a globally coupled system 
is, in general, hard to solve, as the modes may be highly 
singular and it might not be obvious which basis to use 
to expand the functional equation. In the context of the 
model studied in this Letter, we face the additional diffi- 
culty that the microscopic dynamics is characterized by 



a negative Lyapunov exponent and there is no guarantee 
that the distribution P(v,Po,t) is smooth. 
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